#### installing packages if not installed ####

list.of.packages <- c('tidyverse',
                      'ggplot2',
                      'estimatr',
                      'ggthemes',
                      'readr',
                      'readxl',
                      'rdrobust',
                      'gridExtra')
new.packages <- 
  list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### purpose: analyzing philly data ####

#### attaching libraries ####

suppressPackageStartupMessages(
  
  {
    
    library(tidyverse)
    library(ggplot2)
    library(estimatr)    
    library(ggthemes)
    library(readr)
    library(readxl)
    library(rdrobust)
    library(gridExtra)
    
  }
  
)

#### loading data #### 

# stops 

load(file = "data/ph_stops_ped.RData")
load(file = "data/ph_stops_veh.RData")

# crime 

load(file = "data/crimeph.RData")

#### descriptive plots ####

# descriptives: pedestrian stops 

stops_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])
stops_ts = stops_ts %>% filter(blmrv <= 67)

desc_stops1 = stops_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Pedestrian Stops",
       title = "A. Pedestrian Stops") + 
  theme_tufte(base_size = 9)

# descriptives: pedestrian stops 

stops_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])
stops_ts = stops_ts %>% filter(blmrv <= 67)

desc_stops2 = stops_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Vehicle Stops",
       title = "B. Vehicle Stops") + 
  theme_tufte(base_size = 9)

# hit rates (pedestrian)

stops_ped_hit_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            ped_contraband = mean(ped_contraband, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_ped_hit_ts$blmrv = 
  stops_ped_hit_ts$date - min(stops_ped_hit_ts$date[stops_ped_hit_ts$blm == 1])
stops_ped_hit_ts = stops_ped_hit_ts %>% filter(blmrv <= 67)

desc_stops3 = stops_ped_hit_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = ped_contraband),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = ped_contraband, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Hit Rates (Pedestrian)",
       title = "C. Hit Rates (Pedestrian)") + 
  theme_tufte(base_size = 9)

# hit rates (vehicle)

stops_veh_hit_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(veh_contraband, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_veh_hit_ts$blmrv = 
  stops_veh_hit_ts$date - min(stops_veh_hit_ts$date[stops_veh_hit_ts$blm == 1])
stops_veh_hit_ts = stops_veh_hit_ts %>% filter(blmrv <= 67)

desc_stops4 = stops_veh_hit_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = veh_contraband),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = veh_contraband, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Hit Rates (Vehicle)",
       title = "D. Hit Rates (Vehicle)") + 
  theme_tufte(base_size = 9)

# arrest rates (pedestrian)

stops_ped_arr_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            arrest = mean(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_ped_arr_ts$blmrv = 
  stops_ped_arr_ts$date - min(stops_ped_arr_ts$date[stops_ped_arr_ts$blm == 1])
stops_ped_arr_ts = stops_ped_arr_ts %>% filter(blmrv <= 67)

desc_stops5 = stops_ped_arr_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Arrest Rates (Pedestrian)",
       title = "E. Arrest Rates (Pedestrian)") + 
  theme_tufte(base_size = 9)

# arrest rates (vehicle)

stops_veh_arr_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            arrest = mean(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_veh_arr_ts$blmrv = 
  stops_veh_arr_ts$date - min(stops_veh_arr_ts$date[stops_veh_arr_ts$blm == 1])
stops_veh_arr_ts = stops_veh_arr_ts %>% filter(blmrv <= 67)

desc_stops6 = stops_veh_arr_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = arrest),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = arrest, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Arrest Rates (Pedestrian)",
       title = "F. Arrest Rates (Pedestrian)") + 
  theme_tufte(base_size = 9)

# descriptives: black/white stop rate ratio (pedestrian)

stops_ped_rr = stops_ped %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(rr_bw = r_blk_rate / r_wht_rate,
         rr_lw = r_lat_rate / r_wht_rate) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_ped_rr$blmrv = 
  stops_ped_rr$date - min(stops_ped_rr$date[stops_ped_rr$blm == 1])
stops_ped_rr = stops_ped_rr %>% filter(blmrv <= 67)

desc_stops7 = 
  stops_ped_rr %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_bw),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_bw, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "B/W Stop Rate Ratio",
       title = "G. Black/White Stop\nRate Ratio (Pedestrian)") + 
  theme_tufte(base_size = 9)

desc_stops8 = 
  stops_ped_rr %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_lw),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_lw, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "B/L Stop Rate Ratio",
       title = "H. Black/Latino Stop\nRate Ratio (Pedestrian)") + 
  theme_tufte(base_size = 9)

# descriptives: black/white stop rate ratio (vehicle)

stops_veh_rr = stops_veh %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(rr_bw = r_blk_rate / r_wht_rate,
         rr_lw = r_lat_rate / r_wht_rate) %>% 
  filter(date >= as.Date("2020-03-23"))

stops_veh_rr$blmrv = 
  stops_veh_rr$date - min(stops_veh_rr$date[stops_veh_rr$blm == 1])
stops_veh_rr = stops_veh_rr %>% filter(blmrv <= 67)

desc_stops9 = 
  stops_veh_rr %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_bw),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_bw, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "B/W Stop Rate Ratio",
       title = "I. Black/White Stop\nRate Ratio (Vehicle)") + 
  theme_tufte(base_size = 9)

desc_stops10 = 
  stops_veh_rr %>% 
  ggplot() + 
  geom_point(aes(x = date, y = rr_lw),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_lw, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "B/L Stop Rate Ratio",
       title = "J. Black/Latino Stop\nRate Ratio (Vehicle)") + 
  theme_tufte(base_size = 9)

# crime plots 

ph_crime_df_ts = ph_crime_df %>% 
  mutate(date = as.Date(dispatch_date)) %>% 
  dplyr::group_by(date) %>% 
  summarize(crime_violent = sum(crime_violent),
            crime_property = sum(crime_property),
            crime_society = sum(crime_society),
            blm = mean(blm)) %>% 
    filter(date >= as.Date("2020-03-23"))

ph_crime_df_ts$blmrv = 
  ph_crime_df_ts$date - min(ph_crime_df_ts$date[ph_crime_df_ts$blm == 1])
ph_crime_df_ts = ph_crime_df_ts %>% filter(blmrv <= 67)


desc_stops11 = ph_crime_df_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = crime_violent),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_violent, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crime (Against Person)",
       title = "K. Crime\n(Against Person)") + 
  theme_tufte(base_size = 9)

desc_stops12 = ph_crime_df_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = crime_property),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_property, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crime (Against Property)",
       title = "L. Crime\n(Against Property)") + 
  theme_tufte(base_size = 9)

desc_stops13 = ph_crime_df_ts %>% 
  ggplot() + 
  geom_point(aes(x = date, y = crime_society),
             size = .3,
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_society, group = blm),
              se = FALSE,
              col = 'black',
              size = .3) + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             size = .3) + 
  labs(x = "Date", y = "Crime (Against Society)",
       title = "M. Crime\n(Against Society)") + 
  theme_tufte(base_size = 9)

library(grid)

desc_grob = arrangeGrob(desc_stops1, desc_stops2, desc_stops3, desc_stops4,
                        desc_stops5, desc_stops6, desc_stops7, desc_stops8,
                        desc_stops9, desc_stops10, desc_stops11, desc_stops12,
                        desc_stops13, 
                        nrow = 4, ncol = 4)

ggsave(plot = desc_grob, filename = "pics/descph.png", width = 8, height = 8)


#### rdit: cct #### 

# polynomial: 0, 1, 2, 3
# kernel: 0, 1, 2, 3

# ped stop dataset

stops_ts1 = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  mutate(count_ped = count)

stops_ts1$blmrv = 
  stops_ts1$date - min(stops_ts1$date[stops_ts1$blm == 1])

# vehicle stop dataset

stops_ts2 = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  mutate(count_veh = count)

stops_ts2$blmrv = 
  stops_ts2$date - min(stops_ts2$date[stops_ts2$blm == 1])

# ped hit rate

stops_ped_hit_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            ped_contraband = mean(ped_contraband, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ped_hit_ts$blmrv = 
  stops_ped_hit_ts$date - min(stops_ped_hit_ts$date[stops_ped_hit_ts$blm == 1])

# vehicle hit rate

stops_veh_hit_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_contraband = mean(veh_contraband, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_veh_hit_ts$blmrv = 
  stops_veh_hit_ts$date - min(stops_veh_hit_ts$date[stops_veh_hit_ts$blm == 1])

# arrest rates (pedestrian)

stops_ped_arr_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            ped_arrest = mean(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ped_arr_ts$blmrv = 
  stops_ped_arr_ts$date - min(stops_ped_arr_ts$date[stops_ped_arr_ts$blm == 1])

# arrest rates (vehicle)

stops_veh_arr_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            veh_arrest = mean(individual_arrested, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_veh_arr_ts$blmrv = 
  stops_veh_arr_ts$date - min(stops_veh_arr_ts$date[stops_veh_arr_ts$blm == 1])

# rate ratio (pedestrian)

stops_ped_rr_ts = stops_ped %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(ped_rr_bw = r_blk_rate / r_wht_rate,
         ped_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(ped_rr_bw = ifelse(ped_rr_bw == Inf | ped_rr_bw == -Inf, NA, ped_rr_bw),
         ped_rr_lw = ifelse(ped_rr_lw == Inf | ped_rr_lw == -Inf, NA, ped_rr_lw))

stops_ped_rr_ts$blmrv = 
  stops_ped_rr_ts$date - min(stops_ped_rr_ts$date[stops_ped_rr_ts$blm == 1])

# rate ratio (vehicle)

stops_veh_rr_ts = stops_veh %>% 
  group_by(date) %>% 
  summarize(r_blk = sum(r_blk),
            r_lat = sum(r_lat),
            r_wht = sum(r_wht),
            blm = mean(blm)) %>% 
  mutate(r_blk_rate = (r_blk / 614254),
         r_wht_rate = (r_wht / 550102),
         r_lat_rate = (r_lat / 238965)) %>% 
  mutate(veh_rr_bw = r_blk_rate / r_wht_rate,
         veh_rr_lw = r_lat_rate / r_wht_rate) %>% 
  mutate(veh_rr_bw = ifelse(veh_rr_bw == Inf | veh_rr_bw == -Inf, NA, veh_rr_bw),
         veh_rr_lw = ifelse(veh_rr_lw == Inf | veh_rr_lw == -Inf, NA, veh_rr_lw))

stops_veh_rr_ts$blmrv = 
  stops_veh_rr_ts$date - min(stops_veh_rr_ts$date[stops_veh_rr_ts$blm == 1])

# crime 

ph_crime_df_ts = ph_crime_df %>% 
  mutate(date = as.Date(dispatch_date)) %>% 
  dplyr::group_by(date) %>% 
  summarize(crime_violent = sum(crime_violent),
            crime_property = sum(crime_property),
            crime_society = sum(crime_society),
            blm = mean(blm)) 

ph_crime_df_ts$blmrv = 
  ph_crime_df_ts$date - min(ph_crime_df_ts$date[ph_crime_df_ts$blm == 1])

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(stops_ts1 %>% as.data.frame,
                stops_ts2 %>% as.data.frame,
                stops_ped_hit_ts %>% as.data.frame,
                stops_veh_hit_ts %>% as.data.frame,
                stops_ped_arr_ts %>% as.data.frame,
                stops_veh_arr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame)

outcomes = 
  c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
    "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
    "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
    "crime_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                     x = as.numeric(datasets[[i]]$blmrv), 
                     p = polys[j],
                     kernel = kernz[k])
      
      kerns_matrix[k, 1] = out$coef[3]
      kerns_matrix[k, 2] = out$se[3]
      kerns_matrix[k, 3] = out$pv[3]
      kerns_matrix[k, 4] = outcomes[i]
      kerns_matrix[k, 5] = polys[j]
      kerns_matrix[k, 6] = kernz[k]
      
    }
    
    polys_list[[j]] = 
      kerns_matrix %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             outcomes = as.character(outcomes),
             poly = as.numeric(as.character(poly)),
             kern = as.character(kern))
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri.,\n0", "Tri.,\n1", "Tri.,\n2", "Tri.,\n3",
                                      "Uni.,\n0", "Uni.,\n1", "Uni.,\n2", "Uni.,\n3",
                                      "Epa.,\n0", "Epa.,\n1", "Epa.,\n2", "Epa.,\n3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
                                          "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
                                          "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
                                          "crime_society"),
                           labels = c("A. Pedestrian Stops", "A. Traffic Stops",
                                      "C. Hit Rate (Ped)", "B. Hit Rate (Veh)",
                                      "E. Arrest Rate (Ped)", "C. Arrest Rate (Veh)",
                                      "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
                                      "D. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
                                      "E. Crime (Violent)", "F. Crime (Property)", "G. Crime (Society)"))) %>% 
  filter(!grepl(x = outcomes, "edestria|Ped|L\\/W"))

cct_estimates %>% 
  ggplot() + 
  geom_point(aes(x = kernpoly, y = est)) +
  geom_errorbar(aes(x = kernpoly,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se),
                size = .2,
                width = 0) + 
  geom_errorbar(aes(x = kernpoly,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se),
                size = .4,
                width = 0) + 
  geom_hline(yintercept = 0) + 
  facet_wrap(~outcomes, scales = "free", ncol = 4) + 
  labs(x = "Kernel, Polynomial",
       y = "RDiT Coefficient (BLM)") + 
  theme_tufte(base_size = 10) + 
  theme(axis.text.x = element_text(size = 3.5))

ggsave(plot = last_plot(), filename = "pics/cctresph.png", width = 8, height = 3.5)

#### rdit: 10-65 day bandwidth ####

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(stops_ts1 %>% as.data.frame,
                stops_ts2 %>% as.data.frame,
                stops_ped_hit_ts %>% as.data.frame,
                stops_veh_hit_ts %>% as.data.frame,
                stops_ped_arr_ts %>% as.data.frame,
                stops_veh_arr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame)

outcomes = 
  c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
    "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
    "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
    "crime_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    # kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    kerns_list = as.list(rep(NA, length(kernz)))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      
      bw_matrix = matrix(NA, ncol = 7, nrow = length(10:100))
      
      for (l in 1:length(10:100)) {
        
        print(paste0("L iteration ", l))
        out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                       x = datasets[[i]]$blmrv, 
                       p = polys[j],
                       kernel = kernz[k],
                       h = (10:100)[l])
        
        bw_matrix[l, 1] = out$coef[3]
        bw_matrix[l, 2] = out$se[3]
        bw_matrix[l, 3] = out$pv[3]
        bw_matrix[l, 4] = outcomes[i]
        bw_matrix[l, 5] = polys[j]
        bw_matrix[l, 6] = kernz[k]
        bw_matrix[l, 7] = (10:100)[l]
        
      }
      
      kerns_list[[k]] = bw_matrix %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern", "bw")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               outcomes = as.character(outcomes),
               poly = as.numeric(as.character(poly)),
               kern = as.character(kern),
               bw = as.numeric(as.character(bw)))
      
    }
    
    polys_list[[j]] = kerns_list
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri., 0", "Tri., 1", "Tri., 2", "Tri., 3",
                                      "Uni., 0", "Uni., 1", "Uni., 2", "Uni., 3",
                                      "Epa., 0", "Epa., 1", "Epa., 2", "Epa., 3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
                                      "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
                                      "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
                                      "crime_society"),
                           labels = c("A. Pedestrian Stops", "B. Traffic Stops",
                                      "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
                                      "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
                                      "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
                                      "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
                                      "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)"))) 

the_outs = c("A. Pedestrian Stops", "B. Traffic Stops",
             "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
             "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
             "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
             "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
             "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)")

the_outs_labs = c("Pedestrian Stops", "Traffic Stops",
                  "Hit Rate (Ped)", "Hit Rate (Veh)",
                  "Arrest Rate (Ped)", "Arrest Rate (Veh)",
                  "B/W Rate Ratio (Ped)", "L/W Rate Ratio (Ped)",
                  "B/W Rate Ratio (Veh)", "L/W Rate Ratio (Veh)",
                  "Crime (Violent)", "Crime (Property)", "Crime (Society)")

the_outs_plot_list = as.list(rep(NA, length(the_outs_labs)))

for (i in 1:length(the_outs_labs)) {
  
  print(paste0("Iteration ", i))
  
  the_outs_plot_list[[i]] = 
    cct_estimates %>%
    filter(outcomes == the_outs[i]) %>% 
    ggplot() +
    geom_point(aes(x = bw, y = est),
               size = .1,
               alpha = .4) +
    geom_errorbar(aes(x = bw, 
                      ymin = est - 1.96 * se,
                      ymax = est + 1.96 * se),
                  width = 0,
                  size = .1) +
    geom_errorbar(aes(x = bw, 
                      ymin = est - 1.645 * se,
                      ymax = est + 1.645 * se),
                  width = 0,
                  size = .2) + 
    geom_hline(yintercept = 0) + 
    facet_wrap(~ kernpoly, scales = "free") + 
    labs(x = "Bandwidth", y = "RDiT Coefficient\n(BLM Protest)",
         title = the_outs_labs[i]) + 
    theme_tufte(base_size = 9)
  
}

ggsave(plot = the_outs_plot_list[[1]], filename = "pics/c2bwpph1.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[2]], filename = "pics/c2bwpph2.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[3]], filename = "pics/c2bwpph3.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[4]], filename = "pics/c2bwpph4.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[5]], filename = "pics/c2bwpph5.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[6]], filename = "pics/c2bwpph6.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[7]], filename = "pics/c2bwpph7.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[8]], filename = "pics/c2bwpph8.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[9]], filename = "pics/c2bwpph9.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[10]], filename = "pics/c2bwpph10.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[11]], filename = "pics/c2bwpph11.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[12]], filename = "pics/c2bwpph12.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[13]], filename = "pics/c2bwpph13.png", width = 8, height = 6.25)

#### rdit: shaving 100 days after the protest #### 

# loop 

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(stops_ts1 %>% as.data.frame,
                stops_ts2 %>% as.data.frame,
                stops_ped_hit_ts %>% as.data.frame,
                stops_veh_hit_ts %>% as.data.frame,
                stops_ped_arr_ts %>% as.data.frame,
                stops_veh_arr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame)

outcomes = 
  c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
    "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
    "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
    "crime_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    # kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    kerns_list = as.list(rep(NA, length(kernz)))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      
      bw_matrix = matrix(NA, ncol = 7, nrow = length(1:100))
      
      for (l in 1:length(1:100)) {
        
        print(paste0("L iteration ", l))
        dftouse = 
          datasets[[i]] %>% 
          filter(blmrv <= -1 | blmrv >= l)
        
        dftouse$blmrv = 
          ifelse(dftouse$blmrv >= 0, dftouse$blmrv - l, dftouse$blmrv)
        
        out = rdrobust(y = as.numeric(dftouse[, outcomes[i]]),
                       x = dftouse$blmrv, 
                       p = polys[j],
                       kernel = kernz[k])
        
        bw_matrix[l, 1] = out$coef[3]
        bw_matrix[l, 2] = out$se[3]
        bw_matrix[l, 3] = out$pv[3]
        bw_matrix[l, 4] = outcomes[i]
        bw_matrix[l, 5] = polys[j]
        bw_matrix[l, 6] = kernz[k]
        bw_matrix[l, 7] = (1:100)[l]
        
      }
      
      kerns_list[[k]] = bw_matrix %>% 
        as.data.frame %>% 
        `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern", "days_cut")) %>% 
        mutate(est = as.numeric(as.character(est)),
               se = as.numeric(as.character(se)),
               pv = as.numeric(as.character(pv)),
               outcomes = as.character(outcomes),
               poly = as.numeric(as.character(poly)),
               kern = as.character(kern),
               days_cut = as.numeric(as.character(days_cut)))
      
    }
    
    polys_list[[j]] = kerns_list
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri., 0", "Tri., 1", "Tri., 2", "Tri., 3",
                                      "Uni., 0", "Uni., 1", "Uni., 2", "Uni., 3",
                                      "Epa., 0", "Epa., 1", "Epa., 2", "Epa., 3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
                                      "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
                                      "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
                                      "crime_society"),
                           labels = c("A. Pedestrian Stops", "B. Traffic Stops",
                                      "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
                                      "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
                                      "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
                                      "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
                                      "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)"))) 

the_outs = c("A. Pedestrian Stops", "B. Traffic Stops",
             "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
             "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
             "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
             "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
             "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)")

the_outs_labs = c("Pedestrian Stops", "Traffic Stops",
                  "Hit Rate (Ped)", "Hit Rate (Veh)",
                  "Arrest Rate (Ped)", "Arrest Rate (Veh)",
                  "B/W Rate Ratio (Ped)", "L/W Rate Ratio (Ped)",
                  "B/W Rate Ratio (Veh)", "L/W Rate Ratio (Veh)",
                  "Crime (Violent)", "Crime (Property)", "Crime (Society)")

the_outs_plot_list = as.list(rep(NA, length(the_outs_labs)))

for (i in 1:length(the_outs_labs)) {
  
  print(paste0("Iteration ", i))
  
  the_outs_plot_list[[i]] = 
    cct_estimates %>%
    filter(outcomes == the_outs[i]) %>% 
    ggplot() +
    geom_point(aes(x = days_cut, y = est),
               size = .1,
               alpha = .4) +
    geom_errorbar(aes(x = days_cut, 
                      ymin = est - 1.96 * se,
                      ymax = est + 1.96 * se),
                  width = 0,
                  size = .1) +
    geom_errorbar(aes(x = days_cut, 
                      ymin = est - 1.645 * se,
                      ymax = est + 1.645 * se),
                  width = 0,
                  size = .2) + 
    geom_hline(yintercept = 0) + 
    facet_wrap(~ kernpoly, scales = "free") + 
    labs(x = "Days Cut", y = "RDiT Coefficient\n(BLM Protest)",
         title = the_outs_labs[i]) + 
    theme_tufte(base_size = 9)
  
}

ggsave(plot = the_outs_plot_list[[1]], filename = "pics/dcpph1.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[2]], filename = "pics/dcpph2.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[3]], filename = "pics/dcpph3.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[4]], filename = "pics/dcpph4.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[5]], filename = "pics/dcpph5.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[6]], filename = "pics/dcpph6.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[7]], filename = "pics/dcpph7.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[8]], filename = "pics/dcpph8.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[9]], filename = "pics/dcpph9.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[10]], filename = "pics/dcpph10.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[11]], filename = "pics/dcpph11.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[12]], filename = "pics/dcpph12.png", width = 8, height = 6.25)
ggsave(plot = the_outs_plot_list[[13]], filename = "pics/dcpph13.png", width = 8, height = 6.25)

#### rdit placebo tests #### 

# setting up loop

polys = c(0, 1, 2, 3)
kernz = c("uni")
datasets = list(stops_ts1 %>% as.data.frame,
               stops_ts2 %>% as.data.frame,
               stops_ped_hit_ts %>% as.data.frame,
               stops_veh_hit_ts %>% as.data.frame,
               stops_ped_arr_ts %>% as.data.frame,
               stops_veh_arr_ts %>% as.data.frame,
               stops_ped_rr_ts %>% as.data.frame,
               stops_ped_rr_ts %>% as.data.frame,
               stops_veh_rr_ts %>% as.data.frame,
               stops_veh_rr_ts %>% as.data.frame,
               ph_crime_df_ts %>% as.data.frame,
               ph_crime_df_ts %>% as.data.frame,
               ph_crime_df_ts %>% as.data.frame)

outcomes =
 c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
   "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
   "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
   "crime_society")

dataset_list = as.list(rep(NA, length(datasets)))

the_lag_max = as.numeric(as.Date("2020-05-29") - min(datasets[[1]]$date)) - 30

# loop

for (i in 1:length(datasets)) {

 print(paste0("I iteration ", i))
 polys_list = as.list(rep(NA, length(polys)))

 for (j in 1:length(polys)) {

   print(paste0("J iteration ", j))

   kerns_list = as.list(rep(NA, length(kernz)))

   for (k in 1:length(kernz)) {

     print(paste0("K iteration ", j))
     the_lag_max = (as.numeric(as.Date("2020-05-29") -
                                 min(datasets[[i]]$date))) - 30
     the_lag_max2 = as.Date(min(datasets[[i]]$date) + 30)
     lagmat = matrix(NA, ncol = 9, nrow = the_lag_max) %>%
       as.data.frame %>%
       `colnames<-` (c("est", "se", "pv", "outcomes", "poly",
                       "kern", "lag", "lagmax", "lagmaxdate"))

     for (z in 30:the_lag_max) {

       print(paste0("Z iteration ", z))

       datasets[[i]]$blmrv_fake = datasets[[i]]$date -
         datasets[[i]]$date[lead(datasets[[i]]$blmrv, z) == 0 &
                              !is.na(lead(datasets[[i]]$blmrv, z))]

       out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                      x = datasets[[i]]$blmrv_fake,
                      p = polys[j],
                      kernel = kernz[k])

       lagmat[z, 1] = out$coef[3]
       lagmat[z, 2] = out$se[3]
       lagmat[z, 3] = out$pv[3]
       lagmat[z, 4] = outcomes[i]
       lagmat[z, 5] = polys[j]
       lagmat[z, 6] = kernz[k]
       lagmat[z, 7] = z
       lagmat[z, 8] = the_lag_max
       lagmat$lagmaxdate[z] = as.Date(the_lag_max2)

     }

     kerns_list[[k]] = lagmat %>% filter(!is.na(est)) %>%
       mutate(lagmaxdate = as.Date(lagmaxdate, origin = "1970-01-01"))

   }

   polys_list[[j]] =
     kerns_list

 }

 dataset_list[[i]] = polys_list

}

plc_test_out_ph = dataset_list
save(x = plc_test_out_ph, file = "data/plctest_ph.RData")

load("data/plctest_ph.RData")

polys = c(0, 1, 2, 3)
kernz = c("tri", "uni", "epa")
datasets = list(stops_ts1 %>% as.data.frame,
                stops_ts2 %>% as.data.frame,
                stops_ped_hit_ts %>% as.data.frame,
                stops_veh_hit_ts %>% as.data.frame,
                stops_ped_arr_ts %>% as.data.frame,
                stops_veh_arr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame)

outcomes = 
  c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
    "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
    "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
    "crime_society")

dataset_list = as.list(rep(NA, length(datasets)))

for (i in 1:length(datasets)) {
  
  print(paste0("I iteration ", i))
  polys_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J iteration ", j))
    
    kerns_matrix = matrix(NA, ncol = 6, nrow = length(kernz))
    
    for (k in 1:length(kernz)) {
      
      print(paste0("K iteration ", k))
      out = rdrobust(y = as.numeric(datasets[[i]][, outcomes[i]]),
                     x = as.numeric(datasets[[i]]$blmrv), 
                     p = polys[j],
                     kernel = kernz[k])
      
      kerns_matrix[k, 1] = out$coef[3]
      kerns_matrix[k, 2] = out$se[3]
      kerns_matrix[k, 3] = out$pv[3]
      kerns_matrix[k, 4] = outcomes[i]
      kerns_matrix[k, 5] = polys[j]
      kerns_matrix[k, 6] = kernz[k]
      
    }
    
    polys_list[[j]] = 
      kerns_matrix %>% 
      as.data.frame %>% 
      `colnames<-` (c("est", "se", "pv", "outcomes", "poly", "kern")) %>% 
      mutate(est = as.numeric(as.character(est)),
             se = as.numeric(as.character(se)),
             pv = as.numeric(as.character(pv)),
             outcomes = as.character(outcomes),
             poly = as.numeric(as.character(poly)),
             kern = as.character(kern))
    
  }
  
  dataset_list[[i]] = polys_list
  
}

cct_estimates = 
  dataset_list %>% 
  unlist(recursive = FALSE) %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(kernpoly = paste0(kern, ", ", poly)) %>% 
  mutate(kernpoly = factor(kernpoly, 
                           levels = c("tri, 0", "tri, 1", "tri, 2", "tri, 3",
                                      "uni, 0", "uni, 1", "uni, 2", "uni, 3",
                                      "epa, 0", "epa, 1", "epa, 2", "epa, 3"),
                           labels = c("Tri.,\n0", "Tri.,\n1", "Tri.,\n2", "Tri.,\n3",
                                      "Uni.,\n0", "Uni.,\n1", "Uni.,\n2", "Uni.,\n3",
                                      "Epa.,\n0", "Epa.,\n1", "Epa.,\n2", "Epa.,\n3"))) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
                                      "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
                                      "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
                                      "crime_society"),
                           labels = c("A. Pedestrian Stops", "B. Traffic Stops",
                                      "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
                                      "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
                                      "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
                                      "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
                                      "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)"))) 

# placebo test plots 

plc_test_out_ph_df = bind_rows(
  plc_test_out_ph[[1]] %>% bind_rows(),
  plc_test_out_ph[[2]] %>% bind_rows(),
  plc_test_out_ph[[3]] %>% bind_rows(),
  plc_test_out_ph[[4]] %>% bind_rows(),
  plc_test_out_ph[[5]] %>% bind_rows(),
  plc_test_out_ph[[6]] %>% bind_rows(),
  plc_test_out_ph[[7]] %>% bind_rows(),
  plc_test_out_ph[[8]] %>% bind_rows(),
  plc_test_out_ph[[9]] %>% bind_rows(),
  plc_test_out_ph[[10]] %>% bind_rows(),
  plc_test_out_ph[[11]] %>% bind_rows(),
  plc_test_out_ph[[12]] %>% bind_rows(),
  plc_test_out_ph[[13]] %>% bind_rows()
) %>% 
  mutate(outcomes = factor(outcomes,
                           levels = c("count_ped", "count_veh", "ped_contraband", "veh_contraband",
                                      "ped_arrest", "veh_arrest", "ped_rr_bw", "ped_rr_lw",
                                      "veh_rr_bw", "veh_rr_lw", "crime_violent", 'crime_property',
                                      "crime_society"),
                           labels = c("A. Pedestrian Stops", "B. Traffic Stops",
                                      "C. Hit Rate (Ped)", "D. Hit Rate (Veh)",
                                      "E. Arrest Rate (Ped)", "F. Arrest Rate (Veh)",
                                      "G. B/W Rate Ratio (Ped)", "H. L/W Rate Ratio (Ped)",
                                      "I. B/W Rate Ratio (Veh)", "J. L/W Rate Ratio (Veh)",
                                      "K. Crime (Violent)", "L. Crime (Property)", "M. Crime (Society)")))

plc_test_out_ph_df$outcomes_poly = 
  factor(paste0(plc_test_out_ph_df$outcomes, ", ", plc_test_out_ph_df$poly),
         levels = unique(paste0(plc_test_out_ph_df$outcomes, ", ",
                                plc_test_out_ph_df$poly)))

plc_test_out_ph_df$outcomes_poly = 
  factor(plc_test_out_ph_df$outcomes_poly, 
         levels = unique(plc_test_out_ph_df$outcomes_poly))

cct_estimates = cct_estimates %>% filter(kern == "uni")
cct_estimates$outcomes_poly = unique(plc_test_out_ph_df$outcomes_poly)

# before plotting, figuring out the p-value

outcomes_poly_plc = 
  rep(NA, length(unique(plc_test_out_ph_df$outcomes_poly)))


for (i in 1:length(unique(plc_test_out_ph_df$outcomes_poly))) {
  
  print(paste0("Iteration ", i))
  
  outcomes_poly_plc[i] = 
    mean(abs(cct_estimates$est[cct_estimates$outcomes_poly == 
                                 unique(cct_estimates$outcomes_poly)[i]]) >= 
           abs(plc_test_out_ph_df$est[plc_test_out_ph_df$outcomes_poly == 
                                         unique(plc_test_out_ph_df$outcomes_poly)[i]]))
  
}

plc_test_out_ph_df$outcomes_poly = 
  factor(plc_test_out_ph_df$outcomes_poly,
         levels = unique(plc_test_out_ph_df$outcomes_poly),
         labels = paste0(unique(plc_test_out_ph_df$outcomes_poly), ", ",
                         round(outcomes_poly_plc, 2)))

cct_estimates$outcomes_poly = 
  factor(cct_estimates$outcomes_poly,
         levels = unique(cct_estimates$outcomes_poly),
         labels = paste0(unique(cct_estimates$outcomes_poly), ", ",
                         round(outcomes_poly_plc, 2)))

plc_test_out_ph_df = plc_test_out_ph_df %>% 
  filter(!grepl(x = outcomes, 'edestr|Ped|L\\/W')) %>%
  mutate(outcomes_poly = str_replace(outcomes_poly, "^\\w. ", "")) %>% 
  mutate(outcomes_poly = factor(outcomes_poly))

cct_estimates = cct_estimates %>% 
  filter(!grepl(x = outcomes, 'edestr|Ped|L\\/W')) %>% 
  mutate(outcomes_poly = str_replace(outcomes_poly, "^\\w. ", ""))

plc_test_out_ph_df$outcomes_poly = 
  factor(plc_test_out_ph_df$outcomes_poly,
         levels = unique(plc_test_out_ph_df$outcomes_pol))

cct_estimates$outcomes_poly = 
  factor(cct_estimates$outcomes_poly,
         levels = unique(cct_estimates$outcomes_poly))

plc_test_out_ph_df %>% 
  ggplot() + 
  geom_histogram(aes(x = est)) + 
  facet_wrap(~ outcomes_poly, ncol = 4, scales = "free") +
  geom_vline(data = cct_estimates,
             mapping = aes(xintercept = est)) +
  labs(x = "Placebo Coefficients", y = "Count") +  
  theme_tufte(base_size = 8)

ggsave(plot = last_plot(),
       filename = "pics/temp_plc_test_ph.png",
       width = 8, height = 8)

#### saving datasets for all city analysis #### 

datasets_ph = 
  list(stops_ts1 %>% as.data.frame,
                stops_ts2 %>% as.data.frame,
                stops_ped_hit_ts %>% as.data.frame,
                stops_veh_hit_ts %>% as.data.frame,
                stops_ped_arr_ts %>% as.data.frame,
                stops_veh_arr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame,
                stops_ped_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                stops_veh_rr_ts %>% as.data.frame, 
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame,
                ph_crime_df_ts %>% as.data.frame)

save(x = datasets_ph, file = 'data/datasets_ph.Rdata')
